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Abstract 

Here we develop simple numerical algorithms for both stationary and non-stationary solutions of the time- 
dependent Gross-Pitaevskii (GP) equation describing the properties of Bose-Einstein condensates at ultra low tem- 
peratures. In particular, we consider algorithms involving real- and imaginary-time propagation based on a split-step 
Crank-Nicolson method. In a one-space-variable form of the GP equation we consider the one-dimensional, two- 
dimensional circularly-symmetric, and the three-dimensional spherically-symmetric harmonic-oscillator traps. In the 
two-space-variable form we consider the GP equation in two-dimensional anisotropic and three-dimensional axially- 
symmetric traps. The fully-anisotropic three-dimensional GP equation is also considered. Numerical results for the 
chemical potential and root-mean-square size of stationary states are reported using imaginary-time propagation 
programs for all the cases and compared with previously obtained results. Also presented are numerical results of 
non-stationary oscillation for different trap symmetries using real-time propagation programs. A set of convenient 
working codes developed in Fortran 77 are also provided for all these cases (twelve programs in all). In the case of 
two or three space variables, Fortran 90/95 versions provide some simplification over the Fortran 77 programs, and 
these programs are also included (six programs in all). 
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Program summary (1) 

Title of program: imagtimcld.F 

Title of electronic file: imagtimeld.tar.gz 

Catalogue identifier: 

Program summary URL: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland 

Distribution format: tar.gz 

Computers: PC/Linux, workstation/UNIX 

Maximum Ram Memory: 1 GBytc 

Programming language used: Fortran 77 
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Typical running time: Minutes on a medium PC 
Unusual features: None 

Nature of physical problem: This program is designed to solve the time-dependent Gross-Pitacvskii nonlin- 
ear partial differential equation in one space dimension with a harmonic trap. The Gross-Pitacvskii equation 
describes the properties of a dilute trapped Bose-Einstein condensate. 

Method of solution: The time-dependent Gross-Pitaevskii equation is solved by the split-step Crank- 
Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation 
in imaginary time over small time steps. The method yields the solution of stationary problems. 

Program summary (2) 

Title of program: imagtimccir.F 

Title of electronic file: imagtimecir.tar.gz 

Catalogue identifier: 

Program summary URL: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland 

Distribution format: tar.gz 

Computers: PC/Linux, workstation/UNIX 

Maximum Ram Memory: 1 GByte 

Programming language used: Fortran 77 

Typical running time: Minutes on a medium PC 

Unusual features: None 

Nature of physical problem: This program is designed to solve the time-dependent Gross-Pitaevskii non- 
linear partial differential equation in two space dimensions with a circularly-symmetric trap. The Gross- 
Pitacvskii equation describes the properties of a dilute trapped Bose-Einstein condensate. 

Method of solution: The time-dependent Gross-Pitacvskii equation is solved by the split-step Crank- 
Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation 
in imaginary time over small time steps. The method yields the solution of stationary problems. 

Program summary (3) 

Title of program: imagtimcsph.F 

Title of electronic file: imagtimesph. tar.gz 

Catalogue identifier: 

Program summary URL: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland 

Distribution format: tar.gz 

Computers: PC/Linux, workstation/UNIX 

Maximum Ram Memory: 1 GByte 

Programming language used: Fortran 77 

Typical running time: Minutes on a medium PC 

Unusual features: None 

Nature of physical problem: This program is designed to solve the time-dependent Gross-Pitaevskii non- 
linear partial differential equation in three space dimensions with a spherically-symmetric trap. The Gross- 
Pitaevskii equation describes the properties of a dilute trapped Bose-Einstein condensate. 

Method of solution: The time-dependent Gross-Pitaevskii equation is solved by the split-step Crank- 
Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation 
in imaginary time over small time steps. The method yields the solution of stationary problems. 

Program summary (4) 

Title of program: rcaltimeld.F 

Title of electronic file: realtimcld. tar.gz 

Catalogue identifier: 

Program summary URL: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland 

Distribution format: tar.gz 

Computers: PC/Linux, workstation/UNIX 
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Maximum Ram Memory: 2 GBytc 
Programming language used: Fortran 77 
Typical running time: Minutes on a medium PC 
Unusual features: None 

Nature of physical problem: This program is designed to solve the time-dependent Gross-Pitacvskii nonlin- 
ear partial differential equation in one space dimension with a harmonic trap. The Gross-Pitacvskii equation 
describes the properties of a dilute trapped Bose-Einstein condensate. 

Method of solution: The time-dependent Gross-Pitaevskii equation is solved by the split-step Crank- 
Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation 
in real time over small time steps. The method yields the solution of stationary and non-stationary problems. 

Program summary (5) 

Title of program: realtimecir.F 

Title of electronic file: rcaltimccir.tar.gz 

Catalogue identifier: 

Program summary URL: 

Program obtainable from: CPC Program Library. Queen's University of Belfast, N. Ireland 

Distribution format: tar.gz 

Computers: PC/Linux, workstation/UNIX 

Maximum Ram Memory: 2 GByte 

Programming language used: Fortran 77 

Typical running time: Minutes on a medium PC 

Unusual features: None 

Nature of physical problem: This program is designed to solve the time-dependent Gross-Pitaevskii non- 
linear partial differential equation in two space dimensions with a circularly-symmetric trap. The Gross- 
Pitacvskii equation describes the properties of a dilute trapped Bose-Einstein condensate. 

Method of solution: The time-dependent Gross-Pitacvskii equation is solved by the split-step Crank- 
Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation 
in real time over small time steps. The method yields the solution of stationary and non-stationary problems. 

Program summary (6) 

Title of program: realtimesph.F 

Title of electronic file: realtimesph. tar.gz 

Catalogue identifier: 

Program summary URL: 

Program obtainable from: CPC Program Library. Queen's University of Belfast, N. Ireland 

Distribution format: tar.gz 

Computers: PC/Linux, workstation/UNIX 

Maximum Ram Memory: 2 GByte 

Programming language used: Fortran 77 

Typical running time: Minutes on a medium PC 

Unusual features: None 

Nature of physical problem: This program is designed to solve the time-dependent Gross-Pitaevskii non- 
linear partial differential equation in three space dimensions with a spherically-symmetric trap. The Gross- 
Pitaevskii equation describes the properties of a dilute trapped Bose-Einstein condensate. 

Method of solution: The time-dependent Gross-Pitaevskii equation is solved by the split-step Crank- 
Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation 
in real time over small time steps. The method yields the solution of stationary and non-stationary problems. 

Program summary (7) 

Title of programs: imagtimcaxial.F and imagtimcaxial.f90 
Title of electronic file: imagtimeaxial. tar.gz 
Catalogue identifier: 
Program summary URL: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland 
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Distribution format: tar.gz 

Computers: PC/Linux, workstation/UNIX 

Maximum Ram Memory: 2 GByte 

Programming language used: Fortran 77 and Fortran 90 
Typical running time: Few hours on a medium PC 
Unusual features: None 

Nature of physical problem: This program is designed to solve the time-dependent Gross-Pitaevskii non- 
linear partial differential equation in three space dimensions with an axially-symmetric trap. The Gross- 
Pitaevskii equation describes the properties of a dilute trapped Bose-Einstein condensate. 

Method of solution: The time-dependent Gross-Pitaevskii equation is solved by the split-step Crank- 
Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation 
in imaginary time over small time steps. The method yields the solution of stationary problems. 

Program summary (8) 

Title of program: imagtimc2d.F and imagtime2d.f90 
Title of electronic file: imagtime2d. tar.gz 
Catalogue identifier: 
Program summary URL: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland 

Distribution format: tar.gz 

Computers: PC/Linux, workstation/UNIX 

Maximum Ram Memory: 2 GByte 

Programming language used: Fortran 77 and Fortran 90 
Typical running time: Few hours on a medium PC 
Unusual features: None 

Nature of physical problem: This program is designed to solve the time-dependent Gross-Pitaevskii non- 
linear partial differential equation in two space dimensions with an anisotropic trap. The Gross-Pitaevskii 
equation describes the properties of a dilute trapped Bose-Einstein condensate. 

Method of solution: The time-dependent Gross-Pitaevskii equation is solved by the split-step Crank- 
Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation 
in imaginary time over small time steps. The method yields the solution of stationary problems. 

Program summary (9) 

Title of program: realtimeaxial.F and realtimeaxial.f90 
Title of electronic file: realtimeaxial. tar.gz 
Catalogue identifier: 
Program summary URL: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland 

Distribution format: tar.gz 

Computers: PC/Linux, workstation/UNIX 

Maximum Ram Memory: 4 GByte 

Programming language used: Fortran 77 and Fortran 90 
Typical running time: Hours on a medium PC 
Unusual features: None 

Nature of physical problem: This program is designed to solve the time-dependent Gross-Pitaevskii non- 
linear partial differential equation in three space dimensions with an axially-symmetric trap. The Gross- 
Pitaevskii equation describes the properties of a dilute trapped Bose-Einstein condensate. 

Method of solution: The time-dependent Gross-Pitaevskii equation is solved by the split-step Crank- 
Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation 
in real time over small time steps. The method yields the solution of stationary and non-stationary problems. 

Program summary (10) 

Title of program: realtime2d.F and realtime2d.f90 
Title of electronic file: realtime2d. tar.gz 
Catalogue identifier: 
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Program summary URL: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland 

Distribution format: tar.gz 

Computers: PC/Linux, workstation/UNIX 

Maximum Ram Memory: 4 GBytc 

Programming language used: Fortran 77 and Fortran 90 
Typical running time: Hours on a medium PC 
Unusual features: None 

Nature of physical problem: This program is designed to solve the time-dependent Gross-Pitaevskii non- 
linear partial differential equation in two space dimensions with an anisotropic trap. The Gross-Pitaevskii 
equation describes the properties of a dilute trapped Bose-Einstein condensate. 

Method of solution: The time-dependent Gross-Pitaevskii equation is solved by the split-step Crank- 
Nicolson method by discrctizing in space and time. The discretized equation is then solved by propagation 
in real time over small time steps. The method yields the solution of stationary and non-stationary problems. 

Program summary (11) 

Title of program: imagtime3d.F and imagtime3d.f90 
Title of electronic hie: imagtime3d. tar.gz 
Catalogue identifier: 
Program summary URL: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland 

Distribution format: tar.gz 

Computers: PC/Linux, workstation/UNIX 

Maximum Ram Memory: 4 GByte 

Programming language used: Fortran 77 and Fortran 90 
Typical running time: Few days on a medium PC 
Unusual features: None 

Nature of physical problem: This program is designed to solve the time-dependent Gross-Pitaevskii non- 
linear partial differential equation in three space dimensions with an anisotropic trap. The Gross-Pitaevskii 
equation describes the properties of a dilute trapped Bose-Einstein condensate. 

Method of solution: The time-dependent Gross-Pitaevskii equation is solved by the split-step Crank- 
Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation 
in imaginary time over small time steps. The method yields the solution of stationary problems. 

Program summary (12) 

Title of program: rcaltime3d.F and realtimc3d.f90 
Title of electronic file: realtimc3d. tar.gz 
Catalogue identifier: 
Program summary URL: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland 

Distribution format: tar.gz 

Computers: PC/Linux, workstation/UNIX 

Maximum Ram Memory: 8 GByte 

Programming language used: Fortran 77 and Fortran 90 
Typical running time: Days on a medium PC 
Unusual features: None 

Nature of physical problem: This program is designed to solve the time-dependent Gross-Pitaevskii non- 
linear partial differential equation in three space dimensions with an anisotropic trap. The Gross-Pitaevskii 
equation describes the properties of a dilute trapped Bose-Einstein condensate. 

Method of solution: The time-dependent Gross-Pitaevskii equation is solved by the split-step Crank- 
Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation 
in real time over small time steps. The method yields the solution of stationary and non-stationary problems. 
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1. Introduction 

After a successful experimental detection of Bose-Einstein condensates (BEC) of dilute trapped bosonic 
alkali-metal atoms 7 Li, 23 Na, and 87 Rb [1.2] at ultra-low temperatures, there have been intense theoretical ac- 
tivities in studying properties of the condensate using the time-dependent mean-field Gross-Pitaevskii (GP) 
equation under different trap symmetries. Among many possibilities, the following traps have been used in 
various studies: three-dimensional (3D) spherically-symmetric, axially-symmetric and anisotropic harmonic 
traps, two-dimensional (2D) circularly-symmetric and anisotropic harmonic traps, and one-dimensional (ID) 
harmonic trap. The inter-atomic interaction leads to a nonlinear term in the GP equation, which compli- 
cates its accurate numerical solution, specially for a large nonlinearity. The nonlinearity is large for a fixed 
harmonic trap when either the number of atoms in the condensate or the atomic scattering length is large 
and this is indeed so under many experimental conditions. Special care is needed for the solution of the 
time-dependent GP equation with large nonlinearity and there has been an extensive literature on this topic 
[3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,2^^ 

The time-dependent GP equation is a partial differential equation in space and time variables involving 
first-order time and second-order space derivatives together with a harmonic and a nonlinear potential term, 
and has the structure of a nonlinear Schrodingcr equation with a harmonic trap. One commonly used 
procedure for the solution of the time-dependent GP equation makes use of a discretization of this equation 
in space and time and subsequent integration and time propagation of the discretized equation. From a 
knowledge of the solution of this equation at a specific time, this procedure finds the solution after a small 
time step by solving the discretized equation. A commonly used discretization scheme for the GP equation 
is the semi-implicit Crank-Nicolson discretization scheme [49,50,51] which has certain advantages and will 
be used in this work. 

In the simplest one-space-variable form of the GP equation, the solution algorithm is executed in two 
steps. In the first step, using a known initial solution, an intermediate solution after a small interval of time 
A is found neglecting the harmonic and nonlinear potential terms. The effect of the potential terms is then 
included by a first-order time integration to obtain the final solution after time A. In case of two or three 
spatial variables, the space derivatives are dealt with in two or three steps and the effect of the potential 
terms are included next. As the time evolution is executed in different steps it is called a split-step real-time 
propagation method. This method is equally applicable to stationary ground and excited states as well as 
non-stationary states, although in this paper we do not consider stationary excited states. The virtues of 
the semi-implicit Crank-Nicolson scheme [49,50,51] are that it is unconditionally stable and preserves the 
normalization of the solution under real-time propagation. A simpler and efficient variant of the scheme called 
the split-step imaginary-time propagation method obtained by replacing the time variable by an imaginary 
time is also considered. (The GP equation involves complex variables. However, after replacing the time 
variable by an imaginary time the resultant partial differential equation is real, and hence the imaginary- 
time propagation method involves real variables only. This trick leads to an imaginary-time operator which 
results in exponential decay of all states relative to the ground state and can then be applied to any initial 
trial wave function to compute an approximation to the actual ground state rather accurately. We shall 
use imaginary-time propagation to compute the ground state in this paper.) The split-step imaginary-time 
propagation method involving real variables yields very precise result at low computational cost (CPU time) 
and is very appropriate for the solution of stationary problems involving the ground state. The split-step 
real-time propagation method uses complex quantities and yields less precise results for stationary problems; 
however, they are appropriate for the study of non-equilibrium dynamics in addition to stationary problems 
involving excited states also. 

Most of the previous studies [3,4,5,8,9,12,18,20,21,24,25,33,34,35,39,45,46] on the numerical solution of the 
GP equation are confined to a consideration of stationary states only. Some used specifically the imaginary- 
time propagation method [6,31,44,45]. There are few studies [30,37,38,41] for the numerical solution of the 
time-dependent GP equation using the Crank-Nicolson method [49,50,51]. Other methods for numerical solu- 
tion of the time-dependent GP equation have also appeared in the literature [7,14,16,23,26,27,28,29,32,42,43,47,48]. 
These time-dependent methods can be used for studying non-equilibrium dynamics of the condensate in- 
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volving non-stationary states. 

The purpose of the present paper is to develop a simple and efficient algorithm for the numerical solution 
of the GP equation using time propagation together with the semi-implicit Crank-Nicolson discretization 
scheme [49,50,51] specially useful to newcomers in this field interested in obtaining a numerical solution 
of the time-dependent GP equation. Easy-to-use Fortran 77 programs for different trap symmetries with 
adequate explanation arc also included. In case of two and three space variables, Fortran 90/95 programs 
are more compact in nature and these programs are also included. We include programs using both real- and 
imaginary-time propagation. For stationary ground states the imaginary-time method has a much quicker 
convergence rate compared to the real-time method and should be used for the calculation of chemical 
potential, energy and root-mean-square (rms) sizes. We calculate the chemical potential and rms sizes of 
the condensate for stationary problems and compare these results with those previously obtained by other 
workers for different trap symmetries. These results can be easily calculated in a decent PC using the Fortran 
programs provided. In addition to the results for the stationary states, the real-time propagation routines 
can also be used to study the non-stationary transitions, as in collapse dynamics [52] and non-equilibrium 
oscillation [30]. 

This paper is organized as follows. In Sec. 2 we present the GP equations with different traps that we 
consider in this paper, e.g., the 3D spherically-symmetric, 2D circularly-symmetric and ID harmonic traps 
involving one space variable, the anisotropic 2D and axially-symmetric 3D harmonic traps in two space 
variables and the fully anisotropic 3D harmonic trap in three space variables. In Sec. 3 we elaborate the 
numerical algorithm for solving the GP equation in one space variable (the ID, circularly-symmetric 2D, and 
spherically-symmetric 3D cases) and for calculating the chemical potential, energy and rms sizes employing 
both the real- and imaginary-time propagation methods. In Sec. 4 we present the same for solving the GP 
equation in two and three space variables (the anisotropic 2D and axially-symmetric and anisotropic 3D 
cases). In Sec. 5 we present a description of the Fortran programs, an explanation about how to use them, 
and some sample outputs. In Sec. 6 we present the numerical results for chemical potential, rms size, value of 
the wave function at the center, for the ground-state problem using the imaginary-time propagation routines 
and compare our finding with previous results for different trap symmetries in ID, 2D, and 3D. We also 
present a study of non-stationary oscillation in some of these cases using the real-time propagation routines 
when the nonlinear coefficient in the GP equation with a stationary solution was suddenly reduced to half 
its value. Finally, in Sec. 7 we present a brief summary of our study. 



2. Nonlinear Gross-Pitaevskii Equation 

At zero temperature, the time-dependent Bose-Einstein condensate wave function ^ = ^(r; r) at position 
r and time r may be described by the following mean-field nonlinear GP equation [1] 



dfr(r;-r) 
1 dr 



h 2 V 2 

2m 



V(r) + gN\*(r;T)\' 



(1) 



with i = s/—l. Here m is the mass of an atom and N the number of atoms in the condensate, g = Aith 2 a/m 
the strength of inter-atomic interaction, with a the atomic scattering length. The normalization condition 
of the wave function is J c?r|\E'(r; | 2 = 1. 



2.1. Spherically-symmetric GP equation in 3D 

In this case the trap potential is given by V(r) = |rau 2 r 2 , where u> is the angular frequency and f the 
radial distance. After a partial-wave projection the radial part ip of the wave function can be written 
as $(r;r) = ip(f,T). After a transformation of variables to dimensionless quantities defined by r = \plrjl, 
t = Tuj , I = (h/mu>) and (j>(r;t) = <p(r;t)/r = ip(r, r)[/ 3 /(2v / 2)] 1 ^ 2 , the GP equation (1) in this case 
becomes 
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Q2 r 2 



(p{r;t) 



L dt 



¥>(r;t)=0, 



(2) 



where H = %^/2-KNa/l. The purpose of changing the wave function from ip to if = rip is a matter of taste and 
it has certain advantages. First, this transformation removes the first derivative d/dr from the differential 
equation (2) and thus results in a simpler equation [34]. Secondly, at the origin r = 0, ip is a constant, or 
dip/dr = 0. But the new variable satisfies (fi{0,t) = 0. Hence, while solving the differential equation (2), we 
can implement the simple boundary condition that as r — > or oo, ip vanishes. The boundary condition for 
the differential equation in ip will be a mixed one, e.g., the function ip should vanish at infinity and its first 
space derivative should vanish at the origin. The normalization condition for the wave function is 



4tt 



dr\ip(r;t)\ 2 = I. 



(3) 



However, Eq. (2) is not the unique form of dimensionless GP equation in this case. Other forms of dimen- 
sionless equations have been obtained and used by different workers. For example, using the transformations 
r = f/l, t = tuj, I = (h/muj) and (p(r;t) = tp(r;t)/r = ip{f,T)l 3 / 2 , the GP equation (1) becomes 



1 cP_ 

2dr 2 



1 2 

—r 

2 



<p(r;t) 



B 



(4) 



where H = AnNa/l with normalization (3). Finally, using the transformations r 
\J (h/muj) and <p(r;t) = ip(r;t)/r = ip(r,r)l 3 / 2 , the GP equation (1) becomes 



f/l, t = tlj/2, I 



di 



<p(r;t) 



dt 



<p(r;t) = 0, 



(5) 



where K = 8nNa/l with normalization (3). These three sets of dimensionless GP equations have been widely 
used in the literature and will be considered here. Equations (2), (4), and (5) allow stationary solutions 
ip(r;t) = <p(r) exp(— i/j,t) where \i is the chemical potential. The boundary conditions for the solution of 
these equations are ip(0, t) — and linv^oo <p(r, t) = [49]. 



2.2. Anisotropic GP equation in 3D 



The three-dimensional trap potential is given by V(r) = ^mu! 2 {v 2 x 2 + K 2 y 2 + \ 2 z 2 ), where ui x = vuj, lu v = 
luk, and uj z = oj\ are the angular frequencies in the x, y and z directions, respectively, and r = (x, y, z) is the 
radial vector. In terms of dimensionless variables x = \/2x/l, y — \/2y/l, z = \/2z/l, t — tuj, I = \J h/{mui)), 



and (p(x,y,z;t) = J I 3 / (2^/2)^ (r; t), the GP equation (1) becomes 



d 2 d 2 



if- 



dx 2 dy 2 dz 2 ' 4 
with H = 8\/2iraN /I and normalization 



+ - u 2 x 2 + K-y 2 + X 2 z 2 +\t\ip(x,y,z;t)\ 



dt 



<p(x,y,z;t) = o, 



dx 



dy / dz\ip(x,y, z;t)\ 2 = 1. 



(6) 



(7) 



Similarly, using x = x/l,y — y/l,z — z/l,t = tuj, I = y/h/{muj), and tp(x,y, z;t) = Vl 3 ^ (r; r) , the GP 
equation (1) becomes 



Id 2 Id 2 



i d 2 



i 



2dx 2 2dy 2 2 dz 2 2 



- 2 + K 2 y 2 + X 2 z 2 ) +X\^(x,y,z;t)\ 2 -i^ 



tp(x, y, z; t) = 0, (8) 



with H = AiraN/l and normalization (7). Now with scaling t — ► 2t, Eq. (8) can be rewritten as 
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d 2 d 2 d 2 

dx 2 dy 2 dz 2 



+ ( v 2 x 2 + K 2 y 2 + X 2 z 2 I + H \<p(x, y, z; t)\'-i 



ip(x,y,z;t) = 0, 



(9) 



with H = 8naN/l. The boundary conditions for solution are lim-r^^oo (f(x, y, z; t) = 0, lim^-too ip(x, y, z; t) = 
0, lim 2 _ ±00 y>(x, y,z;t) = [49]. 

2.3. Axially- symmetric GP equation in 3D 

In the special case of axial symmetry (y = k) Eqs. (6), (8) and (9) can be simplified considering r = (p, z) 
where p = \/ x 2 + y 2 is the radial coordinate and z is the axial coordinate. Then Eq. (6) becomes 



1 d d 2 



1 



dp 2 p dp dz 2 



+ i "V + AV + h Mp, z- t)r i- 







dt 



(10) 



with H = 8^/2,TTaN/l and normalization 2ir pdp J^ ^ dz\<p(p, z; t)\ 2 = 1. Similarly, Eqs. (8) and (9) can be 
written as 



i d 2 



Id Id 2 



1 



2 dp 2 2pdp 2dz 2 2 
with H = AnaN/l and 

d 2 Id d 2 
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at 



f{p,z;t) = 0, 



9p 2 p 9p <9z 2 



+ H|^,z;t)| 2 -i| 



(p(p,z;t) = 0, 



(11) 



(12) 



with H = 8naN/l. In this case (/?(p = 0, z; t) is not zero but a constant. Convenient boundary conditions for 
solution in this case are lim z ^± 00 <p(p, z\ t) = 0, linip^oo ip(p, z; t) = 0, and dip(p, z; t)/dp\ p —o = [12]. 

2.4. One- dimensional GP equation 

In case of an elongated cigar-shaped trap, which is essentially an axially-symmetric trap with strong trans- 
verse confinement, Eq. (6) reduces to a quasi one-dimensional form. This is achieved by assuming that the 
system remains confined to the ground state in the transverse direction. In this case the wave function of Eq. 
(6) can be written as tp(x, y, z; t) = <p(x; t)(f>o(y)(f)o(z) exp[— i{\+n)t/2] with (j>o(y) = [k/ (27T)] 1 / 4 exp(— ny 2 /4) 
and <po( z ) = [A/ (27T)] 1 / 4 exp(— Az 2 /4) the respective ground state wave functions in y and z directions. Using 
this ansatz in Eq. (6), multiplying by 0o(y)0o( z )j integrating over y and z, dropping the tilde over ip, and 
setting v = 1 we obtain 



d 2 x 2 . . i . . 1 2 . d 



<p{x\t) = 0, 



with H = 2oN\/2\k/1 and normalization 



dx\ip(x; t)\' 



1. 



(13) 



(14) 



Instead if we employ ip(x, y, z; t) = <p(x; i)0o(y)0o(z) exp[— i{\ + n)t/2] with </>o(y) = (n/n) 1 / 4 cxp(— ny 2 /2) 
and .ft, (z) = (A/tt) 1 / 4 cxp(-Az 2 /2) in Eq. (8), in a similar fashion we obtain 



1 t) 2 x 2 1V , , . l2 . d 

-2^ + y +H| ^ ;i)l 



ip(x;t) = 0, 



(15) 



with K = 2aNy\n/l and normalization (14). Now with scaling t — > 2t Eq. (15) can be rewritten as 

-^ + x 2 + K|^;i)| 2 -i- <p(x;t) = 0, (16) 
with H = AaN^/Xn/l and normalization (14). For numerical solution we take lim x _,±oo <p(x,t) = 0. 
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2.5. Anisotropic GP equation in 2D 



In case of a disk-shaped trap, which is essentially an anisotropic trap in two dimensions with strong axial 
binding Eq. (6) reduces to a two-dimensional form. This is achieved by assuming that the system remains 
confined to the ground state in the axial direction. In this case the wave function of Eq. (6) can be written 
as tp(x,y,z;t) — (p(x,y;t)(j)o(z)exp[—iXt/2] with 0o( z ) = [A/(27r)] 1 / 4 cxp(— \z 2 /4) the ground state wave 
function in z direction. Using this ansatz in Eq. (6), multiplying by 0o( z )j integrating over z, dropping the 
tilde over ip and setting v — I wc obtain 



d 2 d 2 x 2 + ny 2 2 . d 



<p(x,y;t) = o, 



now with H = AaN^f2it\/l and normalization 

dx 



dy\tp(x,y]t)\' 



1. 



(17) 



(18) 



Instead if we use in Eq. (8) </>(x, y, z; t) = tf(x,y;t)4>o(z)exp[—iXt/2] with </>q(z) = [X/ir] 1 / 4 cxp(— Az 2 /2), 
then in a similar fashion we obtain 



1 d 2 
Ydx 2 ~ Ydy 2 



1 d 2 x 2 + ny 2 



+ H|^(x,y;t)| 2 -i^ 



f(x,y,t) = 0, 



(19) 



now with N = 2aN^/2i:\/l and normalization (18). Finally, with scaling t — > 2t, Eq. (19) can be written as 

Q2 «2 



^-^+x 2 + K y 2 + ^(x, y ;t)f-i§- t 



v(x,y;t) = o, 



(20) 



with H = 4aN^/2iTX/l and normalization (18). For numerical solution we take \im x ^± oc tp(x, y, t) = and 
]imy^ ±00 tp(x,y,t) = 0. 

2.6. Circularly-symmetric GP equation in 2D 



In the special case of circular symmetry the equations of Sec. 2.5 can be written in one-dimensional form. 
In this case re = 1, and we introduce the radial variable r = (x,y), and rewrite the wave function as ip(r). 
Then the GP equation (17) become 



d 2 Id r 2 , , N l2 d 

The normalization of the wave function is 2tt drr\ip(r; t)\ 2 = 1. 
In the circularly-symmetric case Eq. (19) becomes 

19 2 Id 1. ul . ,. a .8 
W~2^ + r 2 +«l^*>l ~% 



ip(r;t)=0. 



(21) 



ip(r;t)=0. 



Finally, Eq. (20) can be written as 



d 2 Id ^ , . . l2 . d 
^-- 7r +r 2 + H|^(r;t)| 2 -i- 
or z r or at 



ip(r;t) = 0. 



(22) 



(23) 



The convenient boundary condition in this case is lim,.—^ (p(r;t) = and dtp(r;t)/dr\ r= o = [12]. 

In this section we have exhibited GP equations for different trap symmetries. In the next section we 
illustrate the Crank-Nicolson method for the GP equation in one space variable, which is then extended to 
other types of equations in Sec. 4. 
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3. Split-Step Crank-Nicolson Method for the GP Equation in one Space Variable 



3.1. The GP Equation in the ID and radially-symmetric 3D cases 

To introduce the Crank-Nicolson Method [49,50,51] for GP equation we consider first the one-dimensional 
case of Sec. 2.4. The nonlinear GP equation (13) in this case can be expressed in the following form: 



• 9 , , 



d 2 x 2 ... , Sl2 



(p(x;t), 



= H(p(x;t) (24) 

where the Hamiltonian H contains the different linear and nonlinear terms including the spatial derivative. 
(The spherically-symmetric GP equation in 3D has a similar structure and can be treated similarly.) Wc 
solve this equation by time iteration [49,50,51]. A given trial input solution is propagated in time over small 
time steps until a stable final solution is reached. The GP equation is discrctized in space and time using the 
finite difference scheme. This procedure results in a set of algebraic equations which can be solved by time 
iteration using an input solution consistent with the known boundary condition. In the present split-step 
method [50] this iteration is conveniently done in few steps by breaking up the full Hamiltonian into different 
derivative and non-derivative parts. 



3.1.1. Real-time propagation 

The time iteration is performed by splitting H into two parts: H = Hi 

2 



H2, with 



Hi 



H> 



X 
T + 

dx 2 



*\<p{x;t)\' 



Essentially we split Eq. (24) into 



■ 9 . , 



T + «l^i*)l' 



d 

i—^(x-t) = -—^{x;t) 



x 
4 
d 2 
dx 



<p(x; t) = Hiip(x; t) 



H 2 <p(x;t) 



(25) 
(26) 

(27) 
(28) 



Wc first solve Eq. (27) with a initial value <p(x; to) at t = to to obtain an intermediate solution at t — to + A, 
where A is the time step. Then this intermediate solution is used as initial value to solve Eq. (28) yielding 
the final solution at t = to + A as <p(x; to + A). This procedure is repeated n times to get the final solution 
at a given time £fi na i = t,Q + nA. 

The time variable is discretized as t n = nA where A is the time step. The solution is advanced first 
over the time step A at time t n by solving the GP equation (24) with H = Hi to produce an intermediate 
solution (p n+1 / 2 from 93", where <p n is the discretized wave function at time t„. As there is no derivative in 
Hi , this propagation is performed essentially exactly for small A through the operation 



.71+1/2 



<p--- = On d (Hx)<p n = e-^<p n , (29) 

where Ond(-^i) denotes time-evolution operation with Hi and the suffix 'nd' denotes non-derivative. Next 
we perform the time propagation corresponding to the operator H2 numerically by the semi-implicit Crank- 
Nicolson scheme (described below) [49]: 

„n+l _ , n n+l/2 



•r 



-iA 



1+1/2 



The formal solution to (30) is 



n+l 



Ocn(#2)<P 



n+1/2 



l-iAifa/2 
1 + iAi7 2 /2 
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-1/2 



(30) 



(31) 
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which combined with Eq. (29) yields 

V n+1 = OcN(#2)0„d(tfi V\ (32) 
where Ocn denotes time-evolution operation with H2 and the suffix 'CN' refers to the Crank-Nicolson 
algorithm. Operation Ocn 1S used to propagate the intermediate solution tp n+1 / 2 by time step A to generate 
the solution (p n+1 at the next time step t n +i = (n + 1)A. 

The advantage of the above split-step method with small time step A is due to the following three 
factors [50,51]. First, all iterations conserve normalization of the wave function. Second, the error involved 
in splitting the Hamiltonian is proportional to A 2 and can be neglected and the method preserves the 
symplectic structure of the Hamiltonian formulation. Finally, as a major part of the Hamiltonian including 
the nonlinear term is treated fairly accurately without mixing with the delicate Crank-Nicolson propagation, 
the method can deal with an arbitrarily large nonlinear term and lead to stable and accurate converged result. 

Now we describe explicitly the semi-implicit Crank-Nicolson algorithm. The CP equation is mapped onto 
N x one-dimensional spatial grid points in x. Equation (24) is discretized with H = H2 of (26) by the 
following Crank-Nicolson scheme [49,50,51]: 



2h 2 



Grift - 2^ +1 + + fo#i 1/a - 2tf+ 1/a + 



(33) 



where for the spherically-symmetric Eq. (2) = ip(xi;t n ) refers to x = x% = ih, i = 0, 1, 2, N x and h is 
the space step. In the case of the ID Eq. (13), we choose x = Xi = —N x h/2 + ih, i = 0, 1, 2, N x . (The 
choice — even N x — has the advantage of taking an equal number of space points on both sides of x = in 
the ID case setting the point N x /2 at x = 0.) Equation (33) is the explicit form of the formal Eq. (30). This 
scheme is constructed by approximating d/dt by a two-point formula connecting the present (n + 1/2) to 
future (n+1). The spatial partial derivative d 2 /dx 2 is approximated by a three-point formula averaged over 
the present and the future time grid points. This procedure results in a series of tridiagonal sets of equations 
(33) in v?r+i> > an d tPi-i a t time t n +i, which are solved using the proper boundary conditions. 

The Crank-Nicolson scheme (33) possesses certain properties worth mentioning [50,51]. The error in this 
scheme is both second order in space and time steps so that for small A and h the error is negligible. This 
scheme is also unconditionally stable [51]. The boundary condition at infinity is preserved for small values 
of A//i 2 [51]. 

The tridiagonal equations emerging from Eq. (33) are written explicitly as [49] 

AT^+A^ +1 + Attf+{=b u (34) 



where 



A 



— / n+1 12 r, n+1/ 2 . n+l/2x . n+1/ 2 / r\ 

SWi+i +Vi-i ) + <Pi . (35) 



21) 



and A i = Af = — iA/(2/i 2 ), A® = 1 + iA//i 2 . All quantities in hi refer to time step t„+i/2 and are considered 
known. The only unknowns in Eq. (34) are the wave forms <p™±i and </?™ +1 at time step t n +\. To solve Eq. 
(34), we assume the one-term forward recursion relation 

Vi+i = aiLpl +1 +(3 it (36) 

where a, and f3i are coefficients to be determined. Substituting Eq. (36) in Eq. (34) we obtain 

Ar^ + A^l +l + A+( ai tf +1 +0i) = b h (37) 

which leads to the solution 

tf+^JiiArtf+l + A+fr-h), (38) 

with 

7i = -l/(A i+ AtcH). (39) 
From Eqs. (36) and (38) we obtain the following backward recursion relations for the coefficients a, and /3j 

a,_i = 7i Ar, ft^ = 7i (A l +/3 i - h). (40) 
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We shall use the recursion relations (38), (39) and (40) in a backward sweep of the lattice to determine 
on and /3j for i running from N x — 2 down to 0. The initial values chosen are ajv x _i = 0,/?jv x -i = Vjv +1 - 
This ensures the correct value of ip at the last lattice point. After determining the coefficients ai, Pi and 7i, 
we can use the recursion relation (36) from i = to N x — 1 to determine the solution for the entire space 
range using the starting value <^q +1 (=0) known from the boundary conditions. The value at the last lattice 
point is also taken to be known (= 0) . Thus we have determined the solution by using two sets of recursion 
relations across the lattice each involving about N x operations. 

In the numerical implementation of CN real-time propagation the initial state at t = is usually chosen 
to be the analytically known solution of the harmonic potential with zero nonlinearity: H = 0. In the course 
of time iteration the nonlinearity is slowly introduced until the desired final nonlinearity is attained. This 
procedure will lead us to the final solution of the problem. 



3.1.2. Imaginary-time propagation 

Although, real-time propagation as described above has many advantages, in this approach one has to 
deal with complex variables for a complex wave function for non-stationary states. For stationary ground 
state the wave function is essentially real- and the imaginary-time propagation method dealing with real 
variables seems to be convenient. In this approach time t is replaced by an imaginary quantity t = —it and 
Eq. (24) now becomes 

X2 - + K\<p(x;t)f 



dx 2 1 

Hcp(x;t) (41) 



In this equation t is just a mathematical parameter. 

From Eq. (41) we see that an eigenstate ipi of eigenvalue Ei of H satisfying Hipi = Eiipi behaves under 
imaginary-time propagation as d<fi(i)/di = —Ei(pi(i), so that tpi(t) = exp(—Eit)ifi(0). Hence, if we start 
with as arbitrary initial <p(x; t) which can be taken as a linear combination of all cigenfunctions of H, then 
upon imaginary-time propagation all the cigenfunctions will decay exponentially with time. However, all the 
excited states with larger Ei will decay exponentially faster compared to the ground state with the smallest 
eigenvalue. Consequently, after some time only the ground state survives. As all the states are decaying with 
time during imaginary-time propagation, we need to multiply the wave function by a number greater than 
unity to preserve its normalization so that the solution does not go to zero. 

The imaginary-time iteration is performed by splitting H into two parts as before: H = Hi + H2, with Hi 
and H2 given by Eqs. (25) and (26). It is realized that the entire analysis of Sec. 3.1.1 remains valid provided 
we replace i by 1 in Eq. (29) and by —1 in the remaining equations. However, there appears one trouble. 
The CN real-time propagation preserves the normalization of the wave function, whereas the CN imaginary- 
time propagation does not preserve the normalization. This problem can be circumvented by restoring the 
normalization of the wave function after each operation of Crank-Nicolson propagation. Once this is done 
the imaginary-time propagation method for stationary ground state problems yields very accurate result at 
low computational cost. 

Compared to the real-time propagation method, the imaginary-time propagation method is very robust. 
The initial solution in the imaginary-time method could be any reasonable solution and not the analytically 
known solution of a related problem as in the real-time method. Also, the full nonlinearity can be added in 
a small number of time steps or even in a single step and not in a large number of steps as in the real-time 
method. In the programs using the imaginary-time method we include the nonlinearity in a single step. 
These two added features together with the use of real algorithm make the imaginary-time propagation 
method very accurate with quick convergence for stationary ground states as we shall see below. 



3.2. The GP Equation in the circularly-symmetric 2D case 



The Crank-Nicolson discretization for real- and imaginary-time propagation for the circularly-symmetric 
GP equation (21) is performed in a similar fashion as for Eq. (24), apart from the difference that here 
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we also have a first derivative in space variable in addition to the second derivative. Another difference is 
that the wave function is not zero at r = 0. In the case of Eq. (24) we took the boundary condition as 
ip(x;t) = at the boundaries. For the circularly-symmetric case, the convenient boundary conditions are 
lim, — >00 tp(r; t) = and dtp(r; t)/dr\ r= a = 0. 

We describe below the Crank-Nicolson discretization and the solution algorithm in this case for the 
following equation 

_ ld_ _ .d_~ 

dr 2 r dr dt 

required for the solution of Eq. (21). The remaining procedure is similar to that described in detail above 
in the one-dimensional case. 

Equation (42) is discretized by the following Crank-Nicolson scheme as in Eq. (33) [49,50,51]: 



<p(r;t) =0, 



(42) 



n+l 



n+1/2 



1 

2b 2 



w +1 



-1/2 



t+l/2 



l/2^ 



n+1/2 



n+1/2 



(43) 



where again y>™ = 92(1",; t n ), r = r, = ih, ? = 0,1,2,..., N r and h is the space step. This scheme is constructed 
by approximating d/dx by a two-point formula averaged over present and future time grid points. The 
discretization of the first-order time and second-order space derivatives is done as in Eq. (33). This procedure 
results in the tridiagonal sets of equations (43) in (f^+i, an d vf-i at time t n+ i, which are solved 

using the proper boundary conditions. 

The tridiagonal equations emerging from Eq. (43) are written explicitly as 



where 



iA , n+1/2 
^2 



2^ 



-1/2 



1+1/2 



) + </>; 



1+1/2 



iA 

inh 



I n-i 

m+ 



n+1/2 
1 



-1/2, 



(44) 



(45) 



and AT = iA[l/(4/in) - l/(2/i 2 )],A+ = -iA[l / (4h ri ) + l/(2h 2 )}, A° = 1 + iA/h 2 . All quantities in b % refer 
to time step t n +i/2 and are considered known. The only unknowns in Eq. (44) are the wave forms ff+i and 
<p™ +1 at time step t n +\. To solve Eq. (44), we assume the one-term backward recursion relation 



n+l n+l 

fi-i = a m 



■Pi, 



where on and f3i are coefficients to be determined. Substituting Eq. (46) in Eq. (44) we obtain 

A°^ +1 +Ar(a^ +1 +A) = 6 4 , 



which leads to the solution 



with 



bi), 



(46) 



(47) 



(48) 



(49) 



7< = -l/(A i u + A r «i)- 

From Eqs. (46) and (48) we obtain the following forward recursion relations for the coefficients a% and fli 

a i+ i = jiAf, f3 l+1 = jiiArfa - h). (50) 

We shall use the recursion relations (48), (49) and (50) in a forward sweep of the lattice to determine oti and 
Pi for i running from 1 to N r — 1. The initial values chosen are a\ = 1, (3\ = 0. This ensures the correct value 
of the space derivative of ip(r;t) — at r = 0. After determining the coefficients on, Pi and 7.;, we can use 
the recursion relation (46) from i = N r to 1 to determine the solution for the entire space range using the 
starting value ip n+1 (N r ) (=0) known from the boundary condition. Thus we have determined the solution 
by using two sets of recursion relations across the lattice each involving about N r operations. 
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3.3. Chemical potential 



K<p 2 {x) 



(51) 



For stationary states the wave functions for the ID case have the trivial time dependence <p{x; t) 
y>(x) exp(— i/zi), where /i is the chemical potential. Substituting this condition in Eq. (13) we obtain 

x^ 

dx 2 4 

Assuming that the wave form is normalized to unity <p 2 (x)alx = 1, the chemical potential can be 
calculated from the following expression obtained by multiplying Eq. (51) by ip(x) and integrating over all 
space 

2 



/' 



/ dip(a 



\ dx 



+ l p 2 (x){- + ^ 2 (x) 



dx, 



(52) 



where the second derivative has been simplified by an integration by parts. 

All the programs also calculate the many-body energy which is of interest. The analytical expression for 
energy is the same as that of the chemical potential but with the nonlinear term multiplied by 1/2, e.g., [1] 

2 



E 



dip(x) 
dx 



£ 2 (x) 



dx. 



(53) 



The programs will write the value of energy as output. However, we shall not study or tabulate the results for 
energy and we shall not write the explicit algebraic expression of energy in the case of other trap symmetries. 

The GP equation (2) with spherically-symmetric potential is also an one-variable equation quite similar 
in structure to the one-dimensional equation (13) considered above. Hence the entire analysis of Sees. 3.1.1, 
3.1.2, and 3.3 will be applicable in this case with ip(x; t) replaced by ip(r; t)/r in the nonlinear term. Now if we 
consider stationary states of the form <p(r,t) = <p(r) exp(— \\itjYi), the expression for the chemical potential 
becomes 



/i = 4-7T 



/ dipjr) 
V dr 



<p 2 (r) 



dr. 



(54) 



We shall use Eqs. (52) and (54) for the calculation of chemical potential from Eqs. (13) and (2). The energy 
will be calculated from Eq. (53). The expressions for chemical potential in other cases can be written down 
in a straight-forward fashion. 

4. Split-Step Crank-Nicolson method in two and three space variables 

4.1. Anisotropic GP equation in 2D 



In this case the GP equation (17) can be written as 



• 9 , 



= Hip(x,y;t), 



(55) 



with H = 4 v2ttA Na / 1 . The Hamiltonian H can be conveniently broken into three pieces H = H\ + Hi + H3 , 
where 

1 



Hi = - [x 2 + k 2 V 2 J +^{x lV -t)\ 2 , 

H 2 = -—, E* = -—. 

dx 2 ' dy 2 



(56) 
(57) 



Now we adopt a policy quite similar to that elaborated in Sec. 3.1.1 where the Hamiltonian was broken into 
two parts and where the time propagation over time step A using the two parts were carried out alternatively 
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The same procedure will be adopted in the present case where we perform the time propagation using the 
pieces Hi, H2, and H3 of the Hamiltonian successively in three independent time sub-steps A to complete 
a single time evolution over time step A of the entire GP Hamiltonian H. The time propagation over Hi is 
performed as in Eq. (29) and those over H2 and H% as in Eqs. (31) and (33). The chemical potential in this 
case for the stationary state <p(x, y; t) = *p{x, y) exp(— i^t) can be written as 



A* = 



dx 



dy 



dip 



dx ) \ dy 



K 2 1l 2 



(58) 



where we have performed integrations by parts to obtain this form for the chemical potential in terms of 
first derivatives only. 

4.2. Axially- symmetric GP equation in 3D 



In this case the GP equation (10) can be written as 



1 d d 2 



1 



pdp az 2 + l( K V+AV] + %( ft z;t)r 



' dp 2 
= Hip(p,z;t). 

The Hamiltonian H can be conveniently broken into three pieces H = Hi + H2 + H3 , where 

Hi = Up 2 + X 2 z 2 ) + ^(p,z;t)\ 2 , 



(59) 



Ho 



if- 



1 d 



Ho 



(60) 
(61) 



dp 2 pdp' a dz 2 

Now we adopt a policy quite similar to that elaborated in Sec. 4.1 where the Hamiltonian was broken 
into three parts and where the time propagation over time step A using the three parts were carried out 
alternatively. The time propagation over Hi is performed as in Eq. (29) and those over H2 and H3 as in 
Eqs. (43) and (33). The chemical potential in this case for the stationary state <p(p, z; t) = <p(p, z) exp(— ipt) 
can be written as 



pdp 



dz 



dip 



dp 



dip 



dz 



K A p A + \ Z Z 



2 -2 



(62) 



where we have again used integrations by parts to simplify the final expression. 
4.3. Anisotropic GP equation in 3D 



In this case the GP equation (6) can be written as 



• 9 , 



d 2 



d 2 



d 2 



1 •* I 22 ( 22,\22 

~^~7r^^iT2 + -\ vx +K y + A z 

dx dy z dz z 
Hip(x,y,z;t), 



^\ip{x,y,z;t)^ 



ip(x,y,z;t) 



(63) 

with H = 8\/2iraN /I. The Hamiltonian H can be conveniently broken into four pieces H = H1+H2+H3+H4, 
where 



H 1 = ^[v 



Ho 



d 2 



d 2 



Ha 



(64) 
(65) 



d 2 



dx 2 ' dy 2 ' dz 2 

Now we adopt a policy quite similar to that elaborated in Sec. 3.1.1 where the Hamiltonian was broken into 
two parts and where the time propagation over time step A using the two parts were carried out alternatively. 
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The same procedure will be adopted in the present case where we perform the time propagation using the 
pieces Hi , H2 , -H3 and H4 of the Hamiltonian successively in four independent time sub-steps A to complete 
a single time evolution over time step A of the entire GP Hamiltonian H. The time propagation over H\ is 
performed as in Eq. (29) and those over H2, H3 and H4 as in Eqs. (31) and (33). The chemical potential in 
this case for the stationary state ip(x, y, z; t) = tp(x, y, z) exp(— ifit) can be written as 



dx 



dy 



dz 



dip 
dx 



dip 
dy 



d<p 
dz 



v 2 x 2 + K 2 y 2 + X 2 z 2 



Kip' 



(66) 



where we have again used integrations by parts to simplify the final expression. 



5. Description of Numerical Programs 



5.1. GP equation in one space variable 



In this subsection we describe six Fortran codes involving GP equation in one space variable. These are pro- 
grams for solving the ID GP equation (program imagtimeld.F), the circularly-symmetric 2D GP equation 
(program imagtimecir.F) and the radially-symmetric 3D GP equation (program imagtimesph.F) using 
imaginary-time propagation. The similar routines using real-time propagation are realtimeld.F, realtime- 
cir.F and realtimesph.F, respectively. These real- and imaginary-time routines in one space variable have 
similar structures. However, the wave function is real in imaginary-time propagation, whereas it is complex 
in real-time propagation. To accommodate this fact many variables in the real-time propagation routines 
are complex. 

The principal variables employed in the MAIN program arc: N = number of space mesh points, NSTP = 
number of time iterations during which the nonlinearity is introduced in real-time propagation, (in imaginary- 
time propagation the nonlinearity is introduced in one step,) NPAS = number of subsequent time iterations 
with fixed nonlinearity, NRUN = number of final time iterations with (a) fixed nonlinearity in imaginary- 
time propagation and (b) modified nonlinearity to study dynamics in real-time propagation, X(I) = space 
mesh, X2(I) = X(I)*X(I), V(I) = potential, CP(I) = wave function at space point X(I), G, GO = coefficient 
of nonlinear term, MU = chemical potential, EN = energy, ZNORM = normalization of wave function, 
RMS = rms size or radius, DX = space step, DT = time step, OPTION and XOP decide which equation 
to solve. Also important is the subroutine INITIALIZE, where the space mesh X(I), potential V(I), and the 
initial wave function CP (I), are calculated. (An advanced user may need to change the variables V(I) and 
CP (I) so as to adopt the program to solve a different equation with different nonlinearity.) The functions 
and variables not listed above are auxiliary variables, that the user should not need to modify. 

Now we describe the function of the subroutines, which the user should not need to change. The subroutine 
NORM calculates by Simpson's rule the normalization of the wave function and sets the normalization to 
unity, the preassigned value. The real-time propagation preserves the normalization of the wave function 
and hence the subroutine NORM is not used during time propagation. The subroutines RAD and CHEM 
calculate the rms size (length or radius), and chemical potential (and energy) of the system. The subroutines 
COEF and LU together implement the time propagation with the spatial and temporal time derivative 
terms. The subroutine NU performs the time propagation with the nonlinear term and the potential. In the 
imaginary-time program the action of the subroutine LU does not preserve normalization and hence each 
time the subroutine LU is called, the subroutine NORM has to be called to set the normalization of the 
wave function back to unity. (This is not necessary in the real-time programs which preserve the norm.) The 
subroutine NONLIN calculates the nonlinear term. The subroutine DIFF calculates the space derivative 
of the wave function by Richardson's extrapolation formula needed for the computation of the chemical 
potential and energy. The function SIMP does the numerical space integrations with the Simpson's rule. 

The programs implement the splitting method described in Sec. 3.1.1 and 3.2 and calculate the wave 
function, chemical potential, size, normalization, etc. The number of points in the one-dimensional space 
grid represented by the integer variable N has to be chosen consistent with space step DX such that the 
total space covered NxDX is significantly larger than the size of the condensate so that at the boundaries 
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the wave function attains the asymptotic limits (e.g., the absolute value of the wave function or its space 
derivative becomes less than 10 -10 or so for imaginary-time propagation and less than 10 -7 for real-time 
propagation. Note that in the ID and spherically-symmetric 3D problems we are using the asymptotic 
condition that the wave function is zero at both boundaries. In case of the circularly-symmetric 2D problem 
we use a mixed boundary condition, e.g., at origin the derivative of the wave function is zero and at infinity 
the wave function is zero.) In the imaginary-time routine the total space covered should be about 1.5 times 
the extension of the condensate. In the real-time routines, to obtain good precision the total space covered 
should be at least 2 times the extension of the condensate. This is because the imaginary-time routine is 
more precise and the solution attains its limiting asymptotic value at the boundary very rapidly as the total 
space covered is increased. The real-time routine is less accurate and one has to go to a larger distance before 
the solution or its derivative drops to zero. A couple of runs with a sufficiently large N and large DX are 
recommended to have an idea of the size/extension of the system. A smaller value of DX leads to a more 
accurate result provided an appropriate DT is chosen. 

The Crank-Nicolson method, described, for example, by Eq. (33), is unconditionally stable for all A/h 2 . 
Nevertheless, for a numerical application to a specific problem one has to fix the time step DT (= A) and 
space step DX (= h) for good convergence. Space and time steps are given in the MAIN program, as "DATA 
DX/0.0025D0/, DT/0.00002D0/" with correlated DX and DT values obtained by trial. (If the user wants 
to use other values of DX and DT, a set of correlated values obtained by trial is given in Table 1.) The total 
number of space points N in calculation has to be fixed in the line, e.g., "PARAMETER (N = 6000, N2 = 
N/2, NX = N-l)", in the MAIN program, so that the final wave function is within the space range and its 
value is negligibly small at the boundaries. The total numbers of time iterations are fixed in the line, e.g., 
"PARAMETER (NSTP = 500000, NPAS = 10000, NRUN = 20000)", in the MAIN program as described 
below in Sec. 5.3. 

The integer parameter OPTION should be set 1, 2 or 3 in the MAIN routine for solution of equations 
of type (5), (4), or (2), [or for solution of equations of type (16), (15), or (13),] respectively. The difference 
between these three types of equations is in the values of the coefficients in the first two terms. The nonlinear 
term calculated in the subroutine NONLIN is of the form \<p(r, t)\ 2 with coefficient G. If a different type of 
nonlinear term (with different functional dependence on the wave function) is to be introduced, it should 
be done in the subroutine NONLIN. Otherwise, this subroutine should not be changed. A different type of 
nonlinear term is appropriate for a Fermi super-fluid [53,54] or a Tonks-Girardeau gas [55]. 

5.2. GP equation in two and three space variables 

In addition to the programs in one space variable we have six programs in two and three space vari- 
ables. The programs imagtime2d.F and realtime2d.F apply to 2D Cartesian space using imaginary- 
and real-time propagations, respectively. Similarly, imagtime3d.F and realtime3d.F apply to 3D Carte- 
sian space using imaginary- and real-time propagation, respectively. Finally, programs imagtimeaxial.F 
and realtimeaxial.F apply to an axially symmetric trap in 3D. The Fortran 90/95 versions of these pro- 
grams are somewhat more condensed and provide some advantage. Hence we also include these programs 
as imagtime2d.f90, realtime2d.f90, imagtime3d.f90, realtime3d.f90, imagtimeaxial.f90, and real- 
timeaxial.f90. The output of the Fortran 90/95 programs are identical with their corresponding Fortran 
77 versions. All these programs are written using a very similar logic used in the programs in one space 
variable. So most of the considerations described earlier also applies to these cases. We describe the principal 
differences below. 

In case of total number of space points N, one now has the variables NX, NY, and NZ for total number 
of space points in X, Y and Z directions in case of three space variables. In case of two space variables 
the Z component is absent. These variables should be chosen equal to each other for a nearly symmetric 
case. However, if the problem is anisotropic in space, these variables can and should be chosen differently. 
Similarly, the space variable X(I) is now replaced by space variables X(I), Y(J), and Z(K) in three directions. 
Now there are space steps DX, DY, and DZ in place of DX in one space variable. The variables DX, DY, and 
DZ can now be chosen differently in three directions in case of an anisotropic problem. However, they should 
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be chosen together with NX, NY, and NZ so that the wave function lies entirely inside the chosen space and 
becomes negligibly small at the boundaries. If the spatial extension of the wave function is much smaller in 
one space direction than another, one should take a smaller space step in the direction in which the spatial 
extension of the wave function is small. The potential V and wave function CP are now functions of 2 or 
3 space variables and are represented by matrices in place of a column in the case of one space variable. 
The variables AL (or KAP) and BL (or LAM) now define the anisotropy of the trap and are used in the 
subroutine INITIALIZE to define the trap. The subroutine LU is now replaced by LUX, LUY, and LUZ 
to implement the solution in different space directions. In the imaginary-time propagation each time the 
subroutines LUX, LUY, and LUZ operate, one has to set the normalization of the wave function to unity 
by calling the subroutine NORM. 

5.3. Instruction to use the programs 

The programs, as supplied, solve the GP equations for a specific value of nonlincarity and write the wave 
function, chemical potential, energy, rms size or radius, wave function at the center, and nonlincarity, for 
specific values of space and time steps. The real- and imaginary-time programs for a specific equation, for 
example, spherically-symmetric 3D equation, employ similar set of parameters like space and time steps. 
The real-time programs use a larger value of N, so that the discretization covers a larger region in space. 
In all cases the supplied programs solve an equation of type (4) with a factor of 1/2 in front of the space 
derivative, selected by setting the integer parameter OPTION = 2 in the MAIN routine. Other types of 
equations can be obtained by setting OPTION = 1 for equations of type (5) or = 3 for equations of type 
(2). 

For solving a stationary ground state problem, the imaginary-time programs are far more accurate and 
should be used. The real-time programs should be used for studying non-equilibrium problems often using 
an initial wave function calculated by the imaginary-time program. The non-equilibrium problems include 
the study of soliton dynamics [56], expansion [57], collapse dynamics [52], and other types of problems. 

Each program is preset at a fixed nonlinearity GO (= G), correlated DX-DT values and NSTP, NPAS, and 
NRUN. Smaller the steps DX and DT, more accurate will be the result. The correlated values of DX and 
DT on a data line should be found by trial to obtain good convergence. Each supplied program produces 
result up-to a desired precision consistent with the parameters employed — GO, DX, DT, N, NSTP, NPAS, 
and NRUN. If the nonlinearity GO is increased, one might need to increase N to achieve similar precision. 
In many cases one may need an approximate solution (with lower accuracy) involving less CPU time or one 
may need to solve the GP equation for a different value of nonlincarity. 

(a) If GO is reduced, just change the card defining GO. However, if GO is increased, changing the value 
of GO may not be enough. For an increased GO, the wave function extends to a larger region in space. One 
may need to increase the "Number of space mesh points". For a new GO, just plot the output file fort. 3 for 
all programs except realtimc3d.F, realtimc3d.f90, imagtimc3d.F and imagtimc3d.f90, where one should plot 
output files fort.ll, fort. 12, and fort. 13 and see that the wave function is fully and adequately accommodated 
in the space domain. If not, one needs to increase the number in input cards "Number of space mesh points" 
until the wave function is fully and adequately accommodated in the space domain. 

(b) The CPU time involved can be reduced by sacrificing the precision. This can be done by increasing 
the space step(s) and time step and reducing the "Number of space mesh points". If the space step DX 
is increased by a factor of f = 2, the number of space mesh points should be reduced by the same factor. 
The time step DT should be increased by a larger factor, more like f**2. The optimum increase in time 
step should be determined by some experimentation. (A set of correlated DX-DT values is given in Tables 
1 and 4 for the solution of corresponding equations.) The CPU time is the largest in the case of programs 
realtime3d.F, realtime3d.f90, imagtimc3d.F and imagtimc3d.f90 and we give an example of changes below 
in these cases to reduce the CPU time and precision. For example, for realtime3d.F, just change the lines 
15, 16, 17, and 39. The new set of lines should be 

PARAMETER (NX=40, NXX = NX-1, NX2 = NX/2) 
PARAMETER (NY=32, NYY = NY-1, NY2 = NY/2) 
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PARAMETER (NZ=24, NZZ = NZ-1, NZ2 = NZ/2) 
DATA DX /0.5D0/, DY /0.5D0/, DZ /0.5D0/, DT/0.05D0/ 

For imagtimc3d.F, just change the lines 15, 16, 17, and 37. The new set of lines should be 
PARAMETER (NX=24, NXX = NX-1, NX2 = NX/2) 
PARAMETER (NY=20, NYY = NY-1, NY2 = NY/2) 
PARAMETER (NZ=16, NZZ = NZ-1, NZ2 = NZ/2) 
DATA DX /0.5D0/, DY /0.5D0/, DZ /0.5D0/, DT/0.04D0/ 

Please verify, by running the corresponding programs, that the new results so obtained are quite similar 
to the existing results. 

The integer parameter NSTP refers to the number of time iterations during which the nonlinear term is 
slowly introduced in the real-time propagation. This number should be large (typically more than 100,000 
for small nonlincarity, for larger nonlincarity could be 1000,000 ) for good convergence; this means that 
the nonlincarity should be introduced in small amounts over a large number of time iterations. In real-time 
propagation, NPAS refers to certain number of time iterations with the constant nonlinear term already 
introduced in NSTP and should be small (typically 1000). NRUN refers to time iterations with a modified 
nonlincarity so as to generate a non-equilibrium dynamics. In the imaginary-time propagation the parameters 
NPAS and NRUN refer to certain number of time iterations with the constant nonlinear term already 
introduced in one step and should be large (typically NPAS = 200,000 or more) for good convergence. 



5.4. Output files 



The programs write via statements WRITE(1,*), WRITE(2,*) WRITE(3,*) in Files 1, 2, and 3, respec- 
tively, the initial stationary wave function and that after NSTP, and NPAS time iterations for real-time 
propagation and after NPAS and NRUN time iterations for imaginary-time propagation. File 3 gives the 
final stationary wave function of the calculation. However, in the case of the anisotropic 3D programs real- 
time3d.F and imagtime3d.F, sections of the wave function as plotted in Fig. 4 (b) are written in Files 1, 2, 
and 3 before NSTP (realtime3d.F) and after NPAS (imagtime3d.F) iterations, and in Files 11, 12, and 13 
after NPAS (realtime3d.F) and NRUN (imagtime3d.F) iterations. 

In the real-time program a non-stationary oscillation is initiated by suddenly modifying the nonlinearity 
from G to G/2 after NPAS time iterations. During NRUN time iterations the non-stationary dynamics 
is studied. The real-time programs write on File 8 the running time and rms size during non-stationary 
oscillation. 

In addition, these programs write in File 7 the values of nonlincarity G, space steps DX, DY, DZ, time 
step DT, number of space mesh points N, number of time iterations NSTP, NPAS, and NRUN together 
with the values of normalization of the wave function, chemical potential, energy, rms size, value of the wave 
function at center, nonlinearity coefficient G. 

Below wc provide some sample output listed on File 7 from the different programs using OPTION = 2. 
File 7 (fort. 7) represents the comprehensive result in each case. 

(1) Program imagtimeld.F 



OPTION = 2 



# Space Stp N = 8000 

# Time Stp : NPAS = 200000, NRUN = 20000 
Nonlinearity G = 62.74200000 

Space Step DX = 0.002500, Time Step DT = 0.000020 



Norm Chem Ener <x> Psi(0) 

Initial : 1.0000 0.500000 0.500000 0.70711 0.75113 
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After NPAS iter. : 
After NRUN iter . : 



0.9996 
0.9996 



(2) Program imagtimesph.F 
OPTION = 2 



10.369462 
10.369462 



6.256976 
6.256976 



2.04957 
2.04957 



0.40606 
0.40606 



# Space Stp N = 3000 

# Time Stp : , NPAS = 200000, NRUN = 



Nonlinearity G 
Space Step DX = 



20000 



125.48400000 
0.002500, Time Step DT 



0.000020 



Initial : 

After NPAS iter. 

After NRUN iter. 



Norm 



1.0000 
0.9998 
0.9998 



(3) Program imagtimecir.F 
OPTION = 2 



Chem 



Ener 



<r> 



Psi(O) 



1.500000 
4.014113 
4.014113 



1.500000 
3.070781 
3.070781 



1 . 22474 
1.88214 
1.88214 



0.42378 
0. 17382 
0.17382 



# Space Stp N = 2000 

# Time Stp : , NPAS = 200000, NRUN = 200000 



Nonlinearity G 
Space Step DX = 



-2.50970000 
0.002500, Time Step DT = 



0.000020 



Initial : 

After NPAS iter. : 

After NRUN iter. : 



Norm 



1.0000 
1.0000 
1.0000 



Chem 



Ener 



<r> 



Psi(O) 



1.000000 
0.499772 
0.499772 



1.000000 
0.770107 
0.770107 



1.00000 
0.87758 
0.87758 



0.56419 
0.67535 
0.67535 



(4) Program imagtime2d.F and imagtime2d.f90 



OPTION = 2 
Anisotropy AL = 

# Space Stp NX = 

# Time Stp : NPAS 
Nonlinearity G = 
Space Step DX = 
Time Step DT = 



2.000000 

800, NY = 800 

30000, NRUN = 5000 

12.54840000 
0.020000, DY = 0.020000 
0.000100 



Norm 



Chem 



Ener 



<r> 



Psi(0,0) 
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Initial : 1.0000 1.500000 1.500000 0.86603 0.67094 

After NPAS iter. : 0.9999 3.254878 2.490493 1.17972 0.46325 
After NRUN iter.: 0.9999 3.254878 2.490493 1.17972 0.46325 



(5) Program imagtimeaxial.F and imagtimeaxial.fDO 
OPTION = 2 

Anisotropy KAP = 1.000000, LAM = 4.000000 



# Space Stp NX = 500, NY = 500 

# Time Stp : NPAS = 100000, NRUN = 20000 
Nonlinearity G = 18.81000000 

Space Step DX = 0.020000, DY = 0.020000 

Time Step DT = 0.000040 



Norm Chem Energy <rho> <z> <r> psi(O) 



Initial : 1.000 3.00000 3.00000 1.00000 0.35355 1.06066 0.59883 
NPAS iter : 1.000 4.36113 3.78228 1.32490 0.38049 1.37846 0.38129 
NRUN iter : 1.000 4.36113 3.78228 1.32490 0.38049 1.37846 0.38129 



(6) Program imagtime3d.F and imagtime3d.f90 
OPTION = 2 
Anisotropy AL = 1.414214, BL = 2.000000 



# Space Stp NX = 240, NY = 200, NZ = 160 

# Time Stp : NPAS = 5000, NRUN = 500 
Nonlinearity G = 44.90700000 

Space Step DX = 0.050000, DY = 0.050000, DZ = 0.050000 

Time Step DT = 0.000400 



Norm Chem Ener <r> Psi (0,0,0) 



Initial : 1.0000 2.2071 2.2071 1.0505 0.5496 

After NPAS iter. : 0.9997 4.3446 3.4862 1.4583 0.2888 

After NRUN iter. : 0.9997 4.3446 3.4862 1.4583 0.2888 



(7) Program realtimeld.F 
OPTION = 2 

# Space Stp N = 5000 

# Time Stp : NSTP 1000000 , NPAS = 1000, NRUN = 40000 
Nonlinearity G = 62.74200000 

Space Step DX = 0.010000, Time Step DT = 0.000100 



Norm Chem Ener <x> Psi(O) 
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Initial : 1.0000 0.500 0.500 0.707 0.751 

After NSTP iter . : 1.0000 10.368 6.257 2.050 0.406 

After NPAS iter. : 1.0000 10.375 6.257 2.047 0.406 



(8) Program realtimesph.F 
OPTION = 2 

# Space Stp N = 2000 

# Time Stp : NSTP = 1000000, NPAS = 1000, NRUN = 40000 
Nonlinearity G = 125.48400000 

Space Step DX = 0.010000, Time Step DT = 0.000100 



Norm Chem Ener <r> Psi(O) 

Initial : 1.0000 1.500 1.500 1.225 0.424 

After NSTP iter. : 1.0000 4.015 3.071 1.881 0.174 

After NPAS iter. : 1.0000 4.011 3.071 1.884 0.174 



(9) Program realtimecir.F 
OPTION = 2 

# Space Stp N = 2000 

# Time Stp : NSTP = 1000000, NPAS = 1000, NRUN = 40000 
Nonlinearity G = 12.54840000 

Space Step DX = 0.010000, Time Step DT = 0.000100 



Norm 



Chem 



Ener 



<r> 



Psi(O) 



Initial : 

After NSTP iter. 

After NPAS iter. 



1.0000 
1.0000 
1.0000 



1.000 
2.255 
2.255 



1.000 
1.708 
1.708 



(10) Program realtime2d.F and realtime2d.f90 



1.000 
1.308 
1.308 



0.564 
0.391 
0.391 



OPTION = 2 
Anisotropy AL = 

# Space Stp NX = 

# Time Stp : NSTP 
Nonlinearity G = 
Space Step DX = 
Time Step DT = 



1.000000 

200, NY = 200 
100000, NPAS = 1000, NRUN 

12.54840000 
0.100000, DY = 0.100000 
0.001000 



5000 



Norm Chem Ener <r> Psi(0,0) 
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Initial : 1.000 1.000 1.000 1.000 0.564 

After NSTP iter. : 1.000 2.256 1.708 1.307 0.392 

After NPAS iter. : 1.000 2.257 1.708 1.305 0.392 



(11) Program realtimeaxial.F and realtimeaxial.f90 
OPTION = 2 

Anisotropy KAP = 1.000000, LAM = 4.000000 



# Space Stp NX = 130, NY = 130 

# Time Stp : NSTP = 100000, NPAS = 1000, NRUN = 20000 
Nonlinearity G = 18.81000000 

Space Step DX = 0.100000, DY = 0.100000 

Time Step DT = 0.001000 



Norm Chem Energy <rho> <z> psi(0) 



Initial : 1.000 3.000 3.000 1.000 0.354 0.587 

NSTP iter : 0.999 4.362 3.782 1.323 0.381 0.376 

NPAS iter : 1.000 4.362 3.782 1.327 0.379 0.376 



(12) Program realtime3d.F and realtime3d.f90 
OPTION = 2 

Anisotropy AL = 1.414214, BL = 2.000000 



# Space Stp NX = 200, NY = 160, NZ = 120 

# Time Stp : NSTP = 60000, NPAS = 1000, NRUN = 4000 
Nonlinearity G = 22.45400000 

Space Step DX = 0.100000, DY = 0.100000, DZ = 0.100000 

Time Step DT = 0.002000 



Norm Chem Ener <r> Psi (0,0,0) 

Initial : 1.000 2.207 2.207 1.051 0.550 

After NSTP iter. : 1.000 3.572 2.992 1.321 0.347 

After NPAS iter. : 1.000 3.572 2.992 1.320 0.347 
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Fig. 1. (Color online) Relative percentage error as a function of space step for one-dimensional (ID) and spherically-symmetric 
3D (radial) models with nonlincarity K = 627.42 and 627.4, respectively, calculated using the imaginary-time propagation 
programs imagtimeld.F [Eq. (15), OPTION 2] and imagtimesph.F [Eq. (4), OPTION 2] with the data of Table 1. 

6. Numerical Results 

6.1. Stationary Problem 

Table 1 

Convergence of result in the ID (X) and radially-symmetric 3D (r) cases for nonlinearities K = 627.42 and 627.4, calculated 
using the imaginary-time propagation programs imagtimeld.F [Eq. (15), OPTION 2] and imagtimesph.F [Eq. (4), OPTION 
2] respectively, for various space step DX and time step DT. 



H 


DX 


DT 


v(o)Mo) 


^rms / ^rms 


H 


627.42(X) 


0.08 


0.005 


0.276647 


4.384825 


48.024062 


627.42(X) 


0.04 


0.001 


0.276648 


4.384744 


48.024389 


627.42(X) 


0.02 


0.0005 


0.276649 


4.384734 


48.024429 


627.42(X) 


0.01 


0.0001 


0.276649 


4.384726 


48.024462 


627.42(X) 


0.005 


0.00005 


0.276649 


4.384725 


48.024466 


627.42(X) 


0.0025 


0.00002 


0.276649 


4.384724 


48.024468 


627.42(X) 


0.001 


0.00001 


0.276649 


4.384724 


48.024468 


627.4(r) 


0.08 


0.005 


0.106655 


2.506348 


7.247479 


627.4(r) 


0.04 


0.001 


0.106679 


2.505886 


7.248206 


627.4(r) 


0.02 


0.0005 


0.106684 


2.505833 


7.248292 


627.4(r) 


0.01 


0.0001 


0.106686 


2.505785 


7.248365 


627.4(r) 


0.005 


0.00005 


0.106686 


2.505780 


7.248374 


627.4(r) 


0.0025 


0.00002 


0.106686 


2.505776 


7.248380 


627.4(r) 


0.001 


0.00001 


0.106686 


2.505776 


7.248380 



In this subsection we present results for the stationary ground state problem calculated with the imaginary- 
time programs. All numerical results presented in this paper are for the nonlinear equations with a factor 
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Table 2 

The chemical potential fi, rms radius r rms , and the wave function <j>(0) at the center for various nonlincarities in the radially 
symmetric 3D case calculated using the program imagtimcsph.F [Eq. (4), OPTION 2]. Table completed with space step DR 
< 0.0025 and DT = 0.00002. 



N 


0(0) 


0(0) [4] 


''rms 


rrms [4] 




H [3,4] 


-3.1371 


0.48792(1) 


0.4881 


1.51213(1) 


1.1521 


1.265184(2) 


1.2652 





0.42378 


0.4238 


1.22474 


1.2248 


1.500000 


1.5000 


3.1371 


0.38425(1) 


0.3843 


1.27857(1) 


1.2785 


1.677451(1) 


1.6774 


12.5484 


0.31800(1) 


0.3180 


1.39211(1) 


1.3921 


2.065018(1) 


2.0650 


31.371 


0.25810(1) 


0.2581 


1.53561(1) 


1.5356 


2.586116(1) 


2.5861 


125.484 


0.17382(1) 


0.1738 


1.88215(1) 


1.8821 


4.014113(2) 


4.0141 


627.4 


0.10669(1) 


0.1066 


2.50578(1) 


2.5057 


7.248380(3) 


7.2484 


3137.1 


0.06559(1) 


0.0655 


3.41450(1) 


3.4145 


13.553403(4) 


13.553 



of 1/2 in front of the gradient term obtained by choosing OPTION = 2 in the MAIN. First we consider the 
numerical result for the simplest cases — the ID and the radially-symmetric 3D problems by imaginary-time 
propagation with nonlincarities H = 627.42 and 627.4, respectively. The calculations were performed with 
different space and time steps DX and DT, respectively. (As DX is reduced, DT should be reduced also to 
have good convergence. The correlated DX-DT values were obtained by trial to achieve good convergence.) 
In each case a sufficiently large number of space points N is to be taken, so that the space domain of 
integration covers the extension of the wave function adequately. We exhibit in Table 1, the results for the 
wave function at center, rms size, and chemical potential in these two cases for a fixed nonlinearity for 
different DX and DT. We find that convergence is achieved up to six significant digits after the decimal 
point with space step DX = 0.0025 and time step DT = 0.00002. In Fig. 1 we plot the relative percentage 
error in chemical potential \i for various space steps. The percentage error rapidly reduces as space step is 
reduced. 

In Table 2 we exhibit the chemical potential, rms radius, and the wave function at the center for various 
nonlinearities in the spherically symmetric 3D case using space step DX = 0.0025 and time step DT = 0.00002 
in the imaginary-time propagation program. From the table we find that the results are in agreement with 
those calculated in Refs. [3,4]. 

Table 3 

The chemical potential fj,, rms size x rms , and the wave function <p(0) at the center for various nonlinearities in the ID case 
calculated using the program imagtimeld.F [Eq. (15), OPTION 2]. Table completed with DX < 0.0025 and DT = 0.00002. 



N 


tp(0) 


V(0) [4] 


^rms 


Irms [4] 




M [4] 


-2.5097 


0.91317(1) 


0.9132 


0.51334(1) 


0.5133 


-0.80623(3) 


-0.8061 





0.75112 


0.7511 


0.70711 


0.7071 


0.500000 


0.5000 
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0.64596(1) 


0.6459 


0.89602(1) 


0.8960 


1.526593(3) 


1.5265 


12.5484 


0.52975(1) 


0.5297 


1.24549(1) 


1.2454 


3.596560(2) 


3.5965 


31.371 


0.45567(1) 


0.4556 


1.64170(1) 


1.6416 
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6.5526 
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0.4060 
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10.369 
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0.34856(1) 


0.3485 


2.76794(1) 


2.7679 


19.070457(2) 


19.0704 


313.71 


0.31053(1) 


0.3105 


3.48237(1) 


3.4823 


30.259178(3) 


30.259 


627.42 


0.27665(1) 


0.2766 


4.38472(1) 


4.3847 


48.024468(3) 


48.024 


1254.8 


0.24647(1) 


0.2464 


5.52282(1) 


5.5228 


76.226427(3) 


76.226 



In Table 3 we exhibit the chemical potential, rms size, and the wave function at the center for various 
nonlinearities in the ID case using space step DX = 0.0025 and time step DT = 0.00002. From the table we 
find that the results are in agreement with those calculated in Ref. [4]. 
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Fig. 2. (Color online) Plot of wave function profile for the (a) radially-symmetric 3D case [Eq. (4), using imagtimesph.F, 
OPTION 2] and (b) ID case [Eq. (15), using imagtimeld.F, OPTION 2]. The curves are labeled by their respective nonlincaritics 
as tabulated in Tables 2 and 3, respectively. 



Table 4 

Convergence of results for chemical potential fi, rms size r rms and the wave function 0(0) at the center in the Cartesian 2D 
case for a nonlincarity N = 12.5484 and anisotropy k = 1 for different space steps h = DX = DY obtained from the program 
imagtime2d.F [using Eq. (19), OPTION 2]. 





DX=DY 
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^rms 
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1.30678(3) 


2.2559 
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In Figs. 2 (a) and (b) we plot the wave function profiles for the radially-symmetric 3D and ID cases 
calculated using the programs imagtimesph.F and imagtimeld.F, respectively, for different nonlincaritics 
presented in Tables 2 and 3. As the nonlincarity is increased the system becomes more repulsive and the 
wave function extends to a larger domain in space. 

In Table 4 we present the results for the wave function at center, rms radius, and chemical potential 
for the Cartesian 2D case with nonlincarity H = 12.5484 using the imaginary-time propagation program 
imagtime2d.F. We find that desired convergence is achieved with space step DX = 0.02. (Note that the 
converged result in this case is less accurate than those in Tables 1, 2 and 3 as we have used a larger space 
step in order to keep the CPU time small. A finer mesh will increase the accuracy requiring a larger CPU 
time.) 

In Table 5 we present results for the wave function at the center, rms size, and chemical potential for 
different nonlinearities in the anisotropic and circularly-symmetric 2D case calculated using the imaginary- 
time routine imagtimc2d.F and imagtimecir.F, respectively. In the anisotropic case we used space step 0.02 
and time step 0.0001 and in the circularly-symmetric case we used space step 0.0025 and time step 0.00002. 
We also compare these results with those of Rcf. [4] and establish more accurate results in the present 
numerical calculation. 

In Figs. 3 (a) and (b) we plot the wave function profiles for the circularly-symmetric and anisotropic 
2D cases using the programs imagtimecir.F and imagtime2D.F, respectively. In the anisotropic 2D case the 
anisotropy k = 2 and nonlinearity H = 12.4584. Because of anisotropy the wave function in Fig. 3 (b) is 
compressed in the y direction (note the different scales in a; and y directions of the plot.) 
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Table 5 

The chemical potential fj,, rms size r rms , and the wave function ip(0) at the center for various nonlinearities in the anisotropic 
2D case obtained using the programs imagtime2d.F [Eq. (19), OPTION 2] and imagtimecir.F [Eq. (22), OPTION 2], The case 
k = 1 represents circular symmetry and fc ^ 1 corresponds to anisotropy. Table completed with DX = DY < 0.02 and DT 
= 0.0001 for imagtimc2d.F and space step PR < 0.0025 and DT = 0.00002 for imagtimecir.F. 
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1 


627.42 
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Fig. 3. (Color online) Plot of wave function profile for the (a) circularly-symmetric [Eq. (22), OPTION 2] and (b) anisotropic 
2D cases [Eq. (19), OPTION 2] with nonlinearity K = 12.5484 and anisotropy k = 2 using programs imagtimecir.F and 
imagtime2d.F, respectively. Curves in (a) are labeled by the respective nonlinearities as in Table 5. 

In Table 6 we present results for the wave function at the center, rms sizes, and chemical potential 
for different nonlinearities in the axially-symmetric 3D case with n = 1 and A = 4 calculated using the 
imaginary-time routine imagtimeaxial.F and space step Dp = DZ = 0.02 and time step DT = 0.0004. We 
compare with the results of Ref. [4] and establish more accurate results. 

In Figs. 4 (a) and (b) we plot the wave function profiles for the axially-symmetric and fully anisotropic 
3D cases using the programs imagtimeaxial.F and imagtimc3d.F, respectively. In the axially-symmetric 3D 
case the anisotropy k = 1, A = 4 and nonlinearity H = 1881 were employed. In the fully anisotropic 3D case 
the anisotropics v = 1, k = ^2, and A = 2 and nonlinearities H = 22.454, 359.26 and 11496.3 are used. The 
effect of the anisotropy is explicit in the 3D case in generating different profiles of the wave function along 
x, y, and z directions for a fixed nonlinearity. 

Next we consider the fully anisotropic case in three dimensions and compare our results with those of Ref. 
[5] calculated using the program imagtime3d.F. This case mimics a realistic case of experimental interest 
with Na atoms. The completely anisotropic trap considered here is time-orbiting potential (TOP) trap with 
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Table 6 

The chemical potential fi, rms sizes, and the wave function <f>(0) at the center for various nonlinearities in the axially-symmetric 
3D case for 87 Rb atoms in an axial trap with A = 4, k = 1 obtained using the program imagtimeaxial.F [Eq. (11), OPTION 2]. 
The axial frequency is taken as ui z = 80tt Hz, m( 87 Rb) = 1.44 X 10 -25 kg, I = yj h/mui x = 0.3407 X 10~ 5 m, a = 5.1 nm, and 
the ratio between scattering and oscillator length Ana/l = 0.01881. Table completed with DZ < 0.02, Dp < 0.02 and and DT 
= 0.00004. 
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^rms 
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0.5531 
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100000 
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0.1011 


0.1012 


3.2758 


3.276 


0.6173 


0.6174 


19.4751 


19.47 


400000 


7524 


0.0666 


0.0666 


4.3408 


4.341 


0.7881 


0.7881 
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33.47 


800000 
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0.0540 


0.0540 
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0.8976 


0.8976 


44.0234 
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Fig. 4. (Color online) Plot of wave function profile for the (a) axially-symmetric 3D case with nonlincarity K = 1881 anisotropy 
A = 4,k = 1 using imagtimeaxial.F [Eq. (11), OPTION 2] and (b) anisotropic 3D case with nonlincarity K = 22.454,359.26, 
and 11496.3 and anisotropy v = 1,A = \/2 and ft = 2 using imagtime3d.F [Eq. (8), OPTION 2], In the 3D case only the 
sections (f(x, 0, 0), <p(0, y, 0) and <p(0, 0, z) of the wave functions are plotted. 

angular frequencies in the natural ratio (uj x ,ujy,u} z ) = ujq (1, \/2,2), with uj£ = 354tt rad/s. BEC in such 
a system has been observed by Kozuma et al. [58]. We also consider the spherical potential with ojq = 87 
rad/s [59]. The s-wave scattering length of Na is taken as a = 52ao ~ 2.75 nm, with ao = 0.5292 A the Bohr 
radius [5]. 

In Table 7 we exhibit the results for our calculations with the fully anisotropic potential together with 
those for the spherically-symmetric potential in three dimensions. The results for the spherically-symmetric 
potential are also calculated using the one-dimensional radially symmetric imaginary-time program imag- 
timcsph.F in addition to the fully anisotropic program imagtimc3d.F. In the anisotropic case the present 
results are consistent with those of Rcf. [5] . In the spherical case the two sets of the present results (calcu- 
lated with the spherically-symmetric and anisotropic programs) as well as those of Rcf. [5] arc consistent 
with each other. 

The input to our calculation is the nonlincarity coefficient H. which is related to the scattering length a, 
number of atoms N and harmonic oscillator length I via H = AiraN/l in Eq. (8). We provide the nonlincarity 
values of our calculations. Although the present results are in agreement with those of Ref. [5], a very precise 
comparison of the two calculations is not to the point as Schneider and Feder did not provide the nonlincarity 
coefficient H used in their calculation. 
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Table 7 

The chemical potential /i, rms sizes, and wave function </?(0) at the center for various number N of condensate of Na atoms. 
The constants used are m (Na) = 38.175 x 10" 27 kg, a(Na) = 2.75 nm. In all cases the input to numerical calculation was the 
nonlinearity coefficient K = 4irNa/l shown below. For the spherically-symmetric case, we solved the radially symmetric program 
imagtimesph.F [Eq. (4), OPTION 2], (in addition to the 3D anisotropic program imagtime3d.F setting equal frequencies in all 
three directions with DX = DY =DZ = 0.05 and DT = 0.0004) using [5,59] u)j= = 87 rad/s, I = y/h/mu>§ = 5.635 fim, DR. 
< 0.0025 and DT = 0.00002. For the fully anisotropic case we used the program imagtime3d.F [Eq. (8), OPTION 2] with [5,58] 
(j x = u>£ = 354vr rad/s, uj v = V2u x ,u) z =2u> x ,l= yjhfmw§ = 1.576 A*m, DX = DY =DZ = 0.05 and DT = 0.0004. 
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0.2888 


4.3446 


4.345 
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89.81 


1.6328 
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2.970180(1) 


2.9702 
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179.63 


1.8460 
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3.719 
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0.1555 
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4.743445(2) 


4.7434 


4.743 


718.52 


2.3979 
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401.91 


6.123751(2) 


6.1238 


6.124 


1437.03 


2.7447 


0.1022 


15.1284 


15.128 


131072 


803.82 


7.970154(2) 


7.9702 


7.970 


2874.06 


3.1460 


0.0829 


19.8475 


19.847 


262144 


1607.64 


10.426912(3) 


10.4269 


10.427 


5748.13 


3.6092 


0.0673 


26.0961 


26.096 


524288 


3215.28 


13.685486(3) 


13.6855 


13.685 


11496.3 


4.1426 


0.0546 


34.3590 


34.358 



6.2. N on- stationary Oscillation 

The real-time propagation programs calculate the stationary states under different trap symmetries. How- 
ever, they are less efficient than the imaginary-time propagation programs in this task requiring more CPU 
time and producing less accurate results. However, unlike the imaginary-time propagation programs, the 
real-time programs can produce time evolution of non-stationary states also and next we present results of 
such time evolution using the real-time propagation programs under different trap symmetries. 

In this subsection we present results for non-stationary oscillation obtained with the use of the real-time 
programs. After the calculation of the stationary profile, the nonlinearity is suddenly reduced to half. The 
wave function is no longer an cigenstate of the new nonlinear equation. This sets the system into non- 
stationary oscillation which continues for ever. In Fig. 5 we plot the rms size of the wave function vs. time 
t showing this oscillation using the output from File 8 for (a) ID case (using program realtimcld.F), (b) 
radially-symmetric 3D case (using program realtimcsph.F), (c) Cartesian 2D case with anisotropy n = 1 
(using program realtimc2d.F), and (d) 3D axially-symmctric case (using program rcaltimcaxial.F) with 
respective nonlinearitics H = 62.742, 125.484, 12.5484, and 18.81. The rms size at t = is the rms size of 
the stationary wave function obtained after NPAS time iterations. 

Because of transverse instability, the real-time program realtimc3d.F in 3D does not lead to stable si- 
nusoidal oscillation as in other cases for a large change in nonlinearity (nonlinearity reduced to half of its 
initial value) as shown in Fig. 5. Only for small perturbation a sinusoidal oscillation is observed. However, 
we do not present a systematic study of such oscillation. 

7. Summary and Conclusion 

In this paper we describe a split-step method for the numerical solution of the time-dependent nonlinear 
GP equation under the action of a general anisotropic 3D trap using real- and imaginary-time propagation. 
Similar methods for ID and anisotropic 2D traps are also described. The time propagation is carried with 
an initial input. The full Hamiltonian is split into several spatial derivative and a non-derivative parts. The 
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Fig. 5. (Color online) Plot of rms size vs. time for non-stationary oscillation of the system obtained by running the real-time 
programs for (a) ID case (using program realtimeld.F with K = 62.742), (b) spherically-symmetric case (using program 
realtimesph.F with K = 125.484), (c) 2D circularly-symmetric case (using program rcaltimc2d.F with K = 12.5484 and k = 1), 
and (d) 3D axially-symmetric case (using program realtimeaxial.F with K = 18.81 and k = 1, A = 4). The oscillation is started 
during time evolution by suddenly reducing the nonlincarity K to half after the formation of the stationary condensate. 



spatial derivative parts are treated by the Crank-Nicolson method. Different spatial derivative and non- 
derivative parts are dealt in independent steps. This, so called split-step, method leads to highly stable and 
accurate results. 

We considered two types of time iterations — real-time propagation and imaginary-time propagation. In 
the real-time propagation time evolution is performed with the original complex equation. The numerical 
algorithm in this case requires the use of complex variable but produces solution of non-stationary problems. 
In the imaginary-time propagation, the time variable is replaced by i (= s/—T) times a new time variable, 
consequently the GP equation becomes real. The numerical solution of this equation can no longer yield the 
solution of non-stationary problems; but yields very accurate solution of stationary ground state problems 
only, requiring much smaller CPU time. 

We provide the numerical algorithm in detail in ID, 2D, and 3D for real- and imaginary-time propagations. 
We consider six different harmonic oscillator trap symmetries, e.g., a ID trap, a circularly-symmetric 2D trap, 
a radially-symmetric 3D trap, an anisotropic trap in 2D, an axially-symmetric 3D trap, and an anisotropic 
trap in 3D. Each of these cases are treated with real- and imaginary-time propagation algorithms resulting in 
twelve different Fortran 77 programs supplied. We use the imaginary-time propagation programs to provide 
results for different stationary properties of the condensate (chemical potential, rms size, etc) in ID, 2D, 
3D, for different nonlincarities H and compare with previously obtained results [4,5]. In addition we study 
a non-stationary oscillation initiated by suddenly altering the nonlincarity to half its initial value on these 
preformed condensates using the real-time propagation programs. In addition six Fortran 90/95 programs 
are supplied in the case of two and three space variables. 

Although the present programs are valid for the standard GP equation with cubic nonlincarity in a 
harmonic potential, they can be easily adopted for other types of bosonic [55] or fermionic equations [53,54] 
with different nonlinearities and under different types of potentials. To change the potential one should 
change the variable V in the subroutine INITIALIZE and the change in the nonlincarity can be performed 
in the subroutine NONLIN. 
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